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Abstract 


We  develop  a  computationally  efficient  iterative  algorithm  for  source  localization  and 
tracking  using  active/passive  arrays  with  uncertainties  in  sensor  locations.  We  suppose  that 
the  available  data  consist  of  time  delay,  or  differential  time  delay,  measurements  of  the  signal 
wavefront  across  the  array.  We  consider  a  general  senario  in  which  the  array  uncertainties 
may  be  correlated  in  time  and  in  space.  The  proposed  algorithm  is  optimal  in  the  sense  that 
it  converges  montonically  to  the  Maximum  Likelihood  (ML)  estimate  of  the  source  trajectory 
parameters.  In  the  case  of  multiple  sources,  the  algorithm  makes  an  essential  use  of  the 
information  available  from  all  sources  to  reduce  the  array  uncertainties  (the  so-called  array 
calibration)  and  thus  to  improve  the  localization  accuracy  of  each  signal  source.  We  also  de¬ 
rive  new  expressions  for  the  log-likelihood  gradient,  the  Hessian,  and  the  Fisher's  information 
matrix,  that  may  be  used  for  efficient  implementation  of  gradient  based  algorithms,  and  for 
assessing  the  mean  square  error  of  the  resulting  ML  parameter  estimates. 
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I.  Introduction 


The  location  of  a  radiating  source  can  be  determined  by  measurement  of  its  signal  wave- 
front  using  an  array  of  spatially  distributed  sensors/receivers.  In  many  situations  of  practical 
interest,  the  sensor  positions  are  not  precisely  known.  As  demonstrated  in  [l]  [2],  uncertainties 
in  the  sensor  positions  may  seriously  degrade  source  localization  accuracy.  In  such  situations, 
it  has  been  suggested  to  use  auxiliary  sources  at  known/unknown  locations  for  the  purpose  of 
array  calibration  [2]-[5].  However,  most  of  these  studies  have  concentrated  on  analyzing  the 
attainable  mean  square  estimation  errors  without  referring  to  the  estimation  algorithm  that 
may  achieve  the  indicated  performance  predictions. 

In  this  report,  we  develop  a  computationally  efficient  iterative  scheme  for  source  localiza¬ 
tion  and  tracking,  based  on  the  Estimate-Maximize  (EM)  algorithm.  We  suppose  that  the 
available  data  consist  of  time  delay,  or  differential  time  delay,  measurements  of  the  signal 
wavefront  across  the  array.  We  consider  a  general  senario  in  which  the  array  uncertainties 
may  be  correlated  in  time  and  in  space.  The  proposed  algorithm  is  optimal  in  the  sense  that  it 
converges  monotonically  to  the  Maximum  Likelihood  (ML)  estimate  of  the  source  trajectory 
parameters.  In  the  case  of  multiple  sources,  the  algorithm  makes  an  essential  use  of  the  infor¬ 
mation  available  from  all  sources  to  reduce  the  array  uncertainties  (the  so-called  calibration 
effect)  and  thus  to  improve  the  localization  accuracy  of  each  signal  source.  We  also  derive  new 
expressions  for  the  log-likelihood  gradient,  the  Hessian,  and  the  Fisher’s  information  matrix, 
that  may  be  used  for  efficient  implementation  of  gradient-based  algorithms,  and  for  assessing 
the  mean  square  error  (m.s.e.)  of  the  resulting  ML  parameter  estimates. 

The  organization  of  the  paper  is  as  follows:  In  section  II  we  formulate  the  problem  and 
indicate  the  difficulties  associated  with  the  direct  ML  approach.  In  section  III  we  apply  the 
EM  algorithm  for  iteratively  solving  the  ML  problem  in  the  single  source  case,  and  in  section 
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IV  we  consider  the  multiple  source  case.  Section  V  is  devoted  to  gradient-based  algorithms 
and  performance  evaluation,  and  finally  in  section  VI  we  summarize  the  study. 
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II.  Problem  Formulation 


Consider  the  problem  of  source  localization  and  tracking  using  an  array  of  M  spatially 
distributed  receivers  as  illustrated  in  Figure  1.  Let  the  radius  vector  p(i;  0)  denote  the  source 
location  at  time  t,  where  9  is  the  vector  unknown  trajectory  parameters. 

Similarly,  let  qm(t)  denote  the  location  of  the  mth  receiver,  which  is  composed  of  a  nom- 

O 

inaily  known  track  qm  (t),  and  a  random  displacement  MO: 

Qfn(0  =tlm  (0  d"  ^m(0  (1) 


Assuming  perfect  propagation  conditions  in  the  medium,  the  travel  time  of  the  signal  from 
the  source  to  the  mtfl  receiver  is: 


rm[t]  0)  =  -  ||p(i;f)-qm(0  II 

=  i||p(«:»)-4.W-M0 


(2) 


where  ||  •  ||  denotes  the  magnitude  (norm)  of  the  bracketed  vector,  and  c  is  the  propagation 
velocity  in  the  medium.  Using  Taylor  series  expansion  about  the  nominal  sensor  locations: 


*m  —  *m  o  *t*  n 

iq» 


dr„ 


Sm  + 


dh 


m|qm 


2  mdq}mo 


-Sm  +  ' 
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Differentiating  (2), 


Differentiating  (4), 


dr„ 


dq  m  c  ||p  -  qm|| 


;I  — 


,  (p  -  qm)3 


p  -  qmll  Up  -  qmll3 


(p  -  qm)(p  -  qm)3 


where  I  is  the  identity  matrix.  Substituting  (4)  and  (5)  into  (3) 


rm  =  ~  ||  P-  qmll  \ 
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We  suppose  that  the  sensor  displacements  are  small  compared  with  the  spacing  between 
sensors  and  with  the  distance  to  the  source.  Therefore,  ignoring  terms  of  the  order  of 

O 

(||Sm||  /  ||p —  qm  ||)2  (that  is,  the  ratio  of  sensor  displacement  to  source  range  squared),  rm 
is  given  to  an  excellent  approximation  by: 

>»(«;«)  =  ?»(«;*)-;  <£(M)M0  P) 

o 

where  rm  is  the  nominal  delay  measured  relative  to  the  nominal  sensor  position: 

=  (*)!!,  (8) 

and  dm  is  the  unit  vector  along  the  line  connecting  the  nominal  sensor  location  with  the 
source  location  (see  Figure  2): 

llp(‘;*)-  <u  (Oil 

The  product  d^<5m  appearing  in  (7)  is,  therefore,  the  projection  of  5m  in  the  direction  spec¬ 
ified  by  dm.  Geometrically,  this  approximation  means  that  the  signal  wavefront  is  essentially 
planar  over  distances  of  the  order  of  ||$m||. 

Concatenating  the  equations  in  (7)  for  m  =  1,2,...  M ,  we  obtain: 

r(t;0)~?(t;0)-i  DT(t;  fik)6(t)  (10) 

c 

where  r(t;  0)  is  the  vector  time  delays  of  the  signal  wavefront  from  the  source  to  the  actual 
receiver  locations: 

r(t;0)  = 


?2  (t;0) 


(11) 
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r  [t;  0 )  is  the  vector  delays  of  the  signal  from  the  source  to  the  nominal  receiver locations: 

h  (*;#) 

h  m 

(t;  9) 

5(t)  is  the  vector  array  uncertainties: 

h{t) 

m  =  .  (i3) 

Sm(1) 

and  D(t\  0)  is  the  matrix: 

dx(t;*)  (0) 

D{t;0)=  d  2(f,0)  (14) 

■v 

(0)  '  dM(t;0) 

where  the  notation  (0)  indicates  that  all  the  other  components  of  the  matrix  are  zero. 

All  the  information  concerning  source  location  and  track  is  contained  in  the  delays,  or 
differential  delays,  of  the  signed  wavefront  to  the  various  receivers.  Therefore,  source  local¬ 
ization  can  he  performed  by  first  measuring  the  propagation  delays,  and  then  translating  to 
the  estimates  of  the  source  trajectory  parameters.  Since  the  array  uncertainties  enter  only  in 
the  transition  from  the  delay  measurements  to  the  track  estimates,  we  shall  concentrate  only 
on  that  phase  of  the  problem. 

Let  the  vector  delay  measurements  at  time  t  =  tn  be  modeled  as: 

Zn  =  Hn  f(^ni  9)  ~f~  V „ 

~  ff„  ?(*„;«)--  HnDT(tn-,9)8{tn)  +  vn  (15) 

c 
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where  vn  represents  the  measurement  errors.  In  the  transition  from  the  first  line  of  (15)  into 
its  second  line  we  have  invoked  Eq  (10). 

The  matrix  Hn  depends  on  the  measurement  scheme  employed.  If  an  active  system  is 
used  (e.g.,  in  radar),  the  observed  data  consists  of  the  absolute  time-of-arrival  (TOA)  of  the 
signal  wavefront  to  the  various  sensors,  in  which  case  Hn  =  I,  the  idertity  matrix.  In  the 
passive  case  (e.g.,  in  sonar),  the  observed  data  consist  of  time-difference-of-arrival  (TDOA) 
of  the  source  signals  to  the  various  sensor  pairs,  in  which  case  the  matrix  Hn  operates  on 
the  vector  TOA’s  to  generate  the  appropriate  TDOA’s.  For  example,  if  one  measures  the 
TDOA’s  [rv(t)  —  TAf(i)]  *  =  1,2,  ...(M  —  1)  using  sensor  M  as  the  reference,  then  Hn  is  the 
(M  -  1)  x  M  matrix: 

1  -1 

\ 

\  (0)  -1 

'  . 

(0)  \ 

x  t 

\  * 

1  -1 

The  dependence  of  Hn  on  the  time  index  n  accounts  for  situations  in  which  we  may  want 

to  use  different  combinations  of  TOA’s  or  TDOA’s  at  different  times  (perhaps  because  of 
unacceptable  SNR  conditions  at  some  receiver  outputs). 

We  suppose  that  the  measurement  noises  vn  n  =  1,2,...  are  statistically  uncorrelated 
zero-mean  and  Gaussian  with  the  covariance  matrix  Rn ,  depending  on  the  available  signal- 
to-noise  ratio  (SNR)  and  on  the  length  of  the  nt>l  observation  interval.  We  note  that  Rn  is 
generally  a  non-diagonal  matrix  because  of  the  statistical  correlation  between  signal  observa¬ 
tions  at  different  receiver  outputs.  A  lower  bound  on  Rn  that  can  be  used  to  characterize  the 
minimum  attainable  mean  square  error  (m.s.e.)  can  be  found  in  [6], 

We  must  now  make  additional  assumptions  concerning  the  array  uncertainties  S(tn )  n  = 
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1.2.. ..  .  Otherwise,  the  measurement  equation  (15)  contains  too  many  unknowns,  and  the 
problem  cannot  be  solved.  We  suppose  that  S(tn)  n  =  1,2, .. .  satisfy  the  following  stochastic 
difference  equation: 

6{tn)  =  Sn^in-i)  +  Un  n=l,2,...  (16) 

where  6 (to)  is  zero- mean  and  Gaussian  with  a  pre-specified  covariance  matrix  Pq,  and  u„  n  = 

1.2.. ..  are  uncorrelated  zero-mean  and  Gaussian  with  the  covariances  Qn.  This  model  allows 
spatial  as  well  as  temporal  correlation  between  the  array  uncertainties.  We  suppose,  without 
any  loss  of  generality,  that  £(to),  {un},  and  {v„}  are  mutually  independent. 

The  estimation  problem  may  now  be  stated  as  follows:  Given  the  observed  data  z  = 
(zn  :  n  =  1,2, . .  .  N},  find  the  best  possible  estimate  of  the  source  trajectory  parameters  0. 
If  one  interprets  “best”  in  the  usual  sense  of  minimizing  the  mean  square  error  (m.s.e.),  the 
maximum  likelihood  (ML)  estimator  has  a  strong  claim  to  optimality  because  of  its  asymptotic 
efficiency  and  lack  of  bias. 

Observing  the  (15)  subject  to  (16)  constitutes  a  stochastic  dynamic  linear  state  equation 
(where  6(tn)  n  =  0,1,2,...  are  the  state  variables),  the  log-likelihood  function,  that  is  the 
logarithm  of  the  probability  density  of  z,  is  readily  available  in  the  literature  (e.g.,  [7]  [8]): 

=  log/z(M) 

=  6  -  \  J2  {log  det  A"  W  +  #*»  1(®)|*?) }  (17) 

*  n=l 

where  6  is  a  constant  independent  of  0,  and 


M»(«)  =  *n-Hn°r  (tn,0)  +  1  HnDr(tn]  l»)5n/rl_1  (18) 

c 

A ■„(*)  =  ^  HnDT(tn-0)Pn/n,1D{tnJ)HZ  +  Rn  (19) 
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where 


*>n/n- 1  =  E${S(tn)/z  1,  ^2,  ...Zn-l}  (20) 

Pn/n-l  =  E${{6n~  6n/n-i)(6n~  Sn/^^/Zi,  -  -  -  Z„_i }  (21) 

where  1?0  {-/zi,  Z2,  ...  zn-i}  denotes  the  conditioned  expectation  given  Zi,  Z2,  ..  .zn-u  com¬ 
puted  with  respect  to  the  parameter  value  0.  By  applying  the  Kalman  filtering  equations, 
8n/n-i  and  Pn/n- 1  3X6  computed  recursively  in  n. 

Now,  the  ML  estimate  Oml  is  the  0  that  maximizes  L%(0): 

Max  L,(0)  — *  Oml  (22) 

0 

This  is  a  complicated  optimization  problem  that  is  very  difficult  to  solve.  Of  course, 
brute  force  approach  can  always  be  used,  evaluating  the  objective  function  on  a  coarse  grid 
to  roughly  locate  the  global  maximum,  and  then  applying  the  Gauss  method  or  the  Newton- 
Raphson  or  some  other  gradient-search  interative  algorithm.  However,  a  grid  search  involves 
the  computation  of  the  Kalman  filtering  equation  for  a  dense  set  of  parameter  values,  and 
the  computation  of  the  gradient  requires  the  differentiation  of  the  Kalman  filtering  equations. 
Therefore,  these  methods  tend  to  be  computationally  tedius  and  time  consuming. 

In  this  study  we  develop  a  computationally  efficient  iterative  algorithm  for  finding  the  ML 
estimate  of  0,  without  the  need  to  solve  the  complicated  optimization  indicated  above.  The 
key  idea  is  the  if  that  array  uncertainties  <f(tn)  **  =  0,1,2,...  were  known  to  the  observer, 
the  actual  receiver  locations  were  known,  and  the  estimation  problem  would  significantly 
simplifies.  Since  the  array  uncertainties  are  unknown,  they  are  iteratively  estimated  within 
the  algorithm,  and  used  to  iteratively  improve  the  estimate  of  the  various  source  trajectory 
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parameters.  Under  certain  regularity  conditions,  the  proposed  algorithm  converges  monoton- 
ically  to  the  desired  ML  solution. 

First,  we  develop  the  algorithm  for  the  single  source  case,  and  then  present  the  extension 
to  the  multiple  source  case. 
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III.  Development  of  the  Algorithm 

In  this  section  we  apply  the  Estimate-Maximize  (EM)  algorithm  [9]  to  solve  the  indicated 
problem.  The  EM  algorithm  is  basically  an  iterative  method  for  finding  ML  parameter 
estimates.  It  works  with  the  notion  of  complete  data,  and  iterates  between  estimating  the 
log-likelihood  of  the  complete  data  using  the  observed  (incomplete)  data  and  the  current 
parameter  estimate  (El-step),  and  maximizing  the  estimated  log-likelihood  function  to  obtain 
the  new  parameter  estimate  (M-step). 

More  specifically,  let  y  denote  the  complete  data,  that  is  related  to  the  observed  (incom¬ 
plete)  data  z  by  some  non-invertible  (many-to-one)  transformation: 

F(y)  =  z  (23) 

.(a 

Let  9  denote  the  current  estimate  of  0  after  l  iterations  of  the  algorithm.  Then,  the 
next  iteration  cycle  is  specified  in  two  steps  as  follows: 

El-Step:  Compute 


Q(0,d{l))  =  E^{Ly(0)/z}  (24) 

M-Step: 

Max  Q(M(0)  —  $(‘+1)  (25) 

9 

where  Ly{9)  =  log/y(y;0)  is  the  log-likelihood  of  the  complete  data,  and  Egi{-/z} 
denotes  the  conditional  expectation  given  z  evaluated  with  respect  to  the  parameter  value  O' . 

The  EM  algorithm  thus  iteratively  uses  the  latest  parameter  estimate  9^  to  compute 
the  expected  value  of  Ly{9).  It  then  maximizes  over  9  to  get  a  better  parameter  estimate. 
On  the  next  iteration  cycle,  it  uses  the  improved  parameter  estimate  to  improve  the 
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expectation  calculation,  and  thereby  to  improve  the  next  parameter  estimate.  It  has  been 
shown  [9]  [10]  that  if  Q[9,  O')  is  continuous  in  6  and  $' ,  the  EM  algorithm  always  converges  to 
a  stationary  point  of  the  observed  log-likelihood  Lz{0),  where  each  iteration  cycle  increases 
the  value  of  Lz{6).  Of  course,  as  in  all  “hill  climbing”  algorithms,  the  convergence  point  may 
not  be  the  global  maximum,  and  thus  several  starting  points  or  an  initial  grid  search  may  be 
needed. 

To  apply  the  EM  algorithm,  we  need  to  specify  the  complete  data  y.  Clearly,  the  choice  of 
complete  data  is  not  unique  (the  transformation  F(-)  relating  y  to  z  can  be  any  non-invertible 
transformation) .  We  want  to  choose  a  complete  data  y  so  that  the  iterative  expectation  and 
maximization  of  the  associated  log-likelihood  function  Ly  (0)  is  computationally  simpler  than 
the  direct  brute-force  maximixation  of  Lz{6).  As  pointed  out  before,  if  the  array  uncertainties 
were  known  to  the  observer,  the  actual  receiver  locations  where  known,  and  we  would  have 
a  simpler  more  convenient  likelihood  function  to  compute  and  to  maximize.  Therefore,  a 
promising  choice  of  complete  data  y  are  the  delay  observations  z  =  {zn  :  n  =  1,2, . . .  N} 
jointly  with  the  array  uncertainties  S  =  {6(tn)  :  n  =  0, 1,2,...  W}  : 


Applying  Baye’s  rule 

fAT>0)  =  fzlA{z/6).fA{6)  (27) 

where  fA{S)  is  the  probability  density  of  6,  and  fZ/A{z/6)  is  the  conditional  probability 
density  of  z  given  6.  Taking  the  logarithm  on  both  sides  of  (27)  and  involking  (15),  the 
log-likelihood  of  y  is  given  by: 

Ly{9)  =  b  +  \ozfZ/A{z/6) 
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(28) 


n=l 

=*  b-'t\\zn-Hn°r  (t„;0)  +  -cHnDT(tn]0)6{tn)  Hi. 

n=l 

where  b  contains  all  the  terms  that  are  independent  of  0,  including  log  /^(£),  and  where  we 
define: 

II  r  ||ht=  ttW~1t  (29) 

Taking  the  conditional  expectation  of  Ly  (0)  at  a  parameter  value  0^  and  following 
straightforward  matrix  manipulations, 

Q(«,SW)  =  Ejm  {£,(«)/*} 

=  <• -  \  £  {ll  *»  -  nw  f  (lni  +  -HnDr(t„; »)«(')( t„)  ni. 

+±Tr  [/^lffni?r(tn;0)p(1)(in)2?(tn;tf)5T]|  (30) 

where  Tr[-]  stands  for  the  trace  of  the  bracketed  matrix,  and  where  $(0(tn)  and  PM(tn)  are 
the  conditional  mean  estimate  of  8(tn)  and  its  error  convariance  computed  at  $  =  0^: 


W&n)  =  fyo  {S(tn)/Z}  (31) 

P(0(«n)  =  Ef  o  {[«(«„)  ~  5W(t„)]  [fi(tB)  -  (32) 

The  contribution  of  the  Tr[>]  term  appearing  in  (30)  is  very  small  since  P(l)(tn)  is  of  the 
order  of  the  second  moment  (convariance)  of  8(tn).  Therefore,  to  a  very  good  approximation, 

Q(0,0{1))  a  6-|f;i|zn-Fnr(t„;0)  +  iF„Dr(t„;0)^)(t„)||^ 

i  n=l  C 

-  (33) 

L  n=l 


14 


where  in  the  transition  from  the  first  line  of  (33)  to  its  second  line  we  have  used  essentially 
the  same  approximation  as  in  (10).  The  vector 

l)(tn;0) 

0) 

$(*ni*) 

are  the  propagation  delays  computed  with  respect  to  the  estimated  receiver  locations: 

fff(<*;»)  =  i||p(<n;«)-qB(tn)ll  (35) 

c 

where  qm  (*n)  is  the  current  estimate  of  the  mfc^  sensor  location  at  t  =  tn,  which  is  composed 
of  the  nominal  location  compensated  by  the  current  estimate  of  the  sensor  displacement: 

q^=4m(*n) +*&(*«)  (36) 

Recall  (33),  to  compute  Q{0 ,  9^)  (E-Step)  we  only  need  to  compute  S^[t„)  n  =  1,2 
Therefore,  in  this  setting,  the  EM  algorithm  assumes  the  form: 

E-Step:  Compute 

W](tn)  =  EdW{6{tn)/z}  n  =  1,2, ...  iV  (37) 

M-Step: 

Min  £)  ||  *n  -  ,0)  |&,— >  0(‘+1)  (38) 

n=l 

0 

The  algorithm  is  illustrated  in  Figure  3.  We  note  that  the  minimization  in  (38)  is  similar 
to  the  ML  problem  of  estimating  source  trajectory  parameters  when  the  array  uncertainties 
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(sensor  locations)  are  precisely  known  to  the  observer.  The  difference  is  that  in  (38)  we  solve 

the  problem  with  respect  to  the  estimated  sensor  locations.  Thus,  the  algorithm  iterates  back 

.m 

and  forth,  using  the  current  trajectory  estimate  9  to  improve  the  estimate  of  the  array 
uncertainties  in  the  E-Step,  and  thus  to  improve  the  next  estimate  of  the  source  trajectory 
parameters.  Under  the  stated  regularity  condition  (that  is,  the  continuity  of  Q(6,0')  which 
is  clearly  satisfied  here),  the  algorithm  converges  monotonically  to  the  ML  estimate  of  9  or, 
at  least,  to  a  stationary  point  of  the  log-likelihood  function  Lx(9),  where  each  iteration  cycle 
increases  the  likelihood  of  the  estimated  parameters. 

The  conditioned  expectations  required  in  (37)  cam  be  carried  out  efficiently  using  the 
Kalman  smoothing  equations.  Details  are  presented  in  Appendix  A  for  the  convenience  of  the 
reader.  We  have  therefore  substituted  the  direct  ML  approach  that  requires  the  computation 
of  the  Kadman  filtering  equations  at  a  dense  set  of  paramter  values,  with  an  iterative  procedure 
that  involves  the  Kalman  smoothing  equations. 

The  computation  of  6^(£n)  simplifies  if  6(tn)  n  =  1,2,**«  aue  statistically  independent. 
This  is  a  special  case  of  the  model  in  (16)  where  =  0,  amd  it  corresponds  to  a  situation  in 
which  the  time  differences  ( t„  —  tn-i)  exceed  the  correlation  time  of  the  uncertadnty  process. 
Recall  (15),  if  6{tn)  n  =  1,2,-  ••  are  statistically  independent  then  the  pair  (zn,  £(£„))  is 
statisticadly  independent  of  zm  m  n.  Therefore,  to  compute  6^\tn),  we  only  need  to 
compute  the  conditional  expectation  of  6(tn)  with  respect  to  z„.  Since  zn  and  6{tn)  are 
jointly  Gaussian,  satisfying  the  linear  relation  in  (15),  this  conditional  expectation  is  easily 
evaduated: 

6{l)(tn)  =  fyl)  {*(*»)/*.} 

=  --QnD(tn;  9{l))H* 
c 


\HnDT(tn-,d{l))QnD(tn;e{t))HZ  +  Bn 
c 
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(39) 


zn-Hnr(tn;0{l)) 

where  Qn  is  the  a-priori  convariance  of  6  (tn) .  Since  6^  ^  (tn)  depends  only  on  the  measurements 
zn  at  t  =  tn,  it  can  be  computed  recursively  in  n.  There  is  no  need  to  store  past  measurements 
for  the  purpose  of  estimating  the  array  uncertainties. 
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IV.  Multiple  Sources 


We  now  extend  the  scope  by  considering  an  array  of  M  spatially  distributed  sensors, 
receiving  signals  from  K  spatially  distributed  sources,  as  illustrated  in  Figure  4.  Let  p(t;  Ok) 
denote  the  trajectory  of  the  kth  signal  source.  More  generally,  we  could  use  different  functions 
Pk(t\0k)  ta  characterize  the  different  tracks. 

Let  the  delay  measurements  associated  with  the  kth  signal  source  at  t  =  tn  be  given  by: 

zW  =  Ok)  +  vW 

«  T  (tn;  Ok)  -  -H^DT(tn-,  0k)6{tn)  +  yW  (40) 

c 

where  r(t;0k)  is  the  vector  propagation  delays  measured  from  the  kth  source  to  the  actual 
reveiver  locations,  r  (t; &t)  is  the  vector  delays  measured  relative  to  the  nominal  receiver 
locations,  D{t;  •)  is  the  matrix  defined  in  (14),  and  i(tn)  is  the  vector  array  uncertainties 
satisfying  the  stochastic  model  in  (16). 

We  suppose  that  the  measurement  noises  Vn^  associated  with  different  sources  are  mutu¬ 
ally  independent  zero-mean  and  Gaussian  with  the  convariance  This  assumption  means 
that  the  signals  observed  from  different  sources  do  not  interfere  with  each  other  (e.g.,  by 
confining  them  to  disjoint  frequency  bands). 

The  estimation  problem  may  be  s  ated  as  follows:  Given  the  observed  data  z^  k  = 
1,2,.. .  Kt  n  =  1,2, ...  N,  find  the  ML  estimate  of  the  vector  parameters: 

(41) 

Since  the  various  z are  statistically  correlated  due  to  the  array  uncertainties,  the  log- 


0 1 
O2 

Ok 
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likelihood  function  is  now  even  more  complicated,  and  its  maximization  with  respect  to  all  the 


components  of  9  is  a  prohibitive  task.  Of  course,  we  could  decompose  this  complicated  multi¬ 
dimensional  optimization  by  processing  each  delay  equation  separately  to  obtain  the  estimate 
of  the  corresponding  source  location  and  track.  However,  by  that  we  completely  ignore  the 
geometrical  constraint  imposed  by  the  various  sources  to  reduce  the  array  uncertainties  (the 
so-called  calibration  effect  [2]  [3]  [5])  and  thus  to  improve  the  localization  accuracy  of  each 
signal  source. 

We  want  to  apply  the  EM  algorithm  to  solve  the  problem.  As  before,  let  the  complete 
data  y  be  composed  of  the  delay  observations  z  =  j  z^  j  jointly  with  the  array  uncertainties 
6  =  {J(tn)}.  Recall  (40),  and  invoking  the  statistical  independence  of  the  various  vj^’s,  the 
log- likelihood  of  y  is  given  by: 


Ly{9)  =  log  /y(y;0) 

=  log  fA{6 )  +  log  f%/A{z/S) 

*  *=1 n=l  " 

~  b  ~  \  E  E  H  W  °T  (*«:'*)  +  lrtk)VT(tn;  h)6(tn)  |£<*,  (42) 

1  n~  1 

where  6  contains  all  the  terms  that  are  independent  of  0. 

Taking  the  conditional  expectation  of  Ly(9)  at  a  parameter  value  and  ignoring  terms 
of  the  order  of  0  (||  6(tn)  ||2)  (see  the  considerations  leading  from  (28)  to  (33)), 


Q(MW)  =  EdW{Ly{0)/z} 


K  S 


b  -  2  E  £  II  ^ '  -  H{nk)  r  (*.;•*)  +  -HWDT(tn]9k)6{l\tn)  ||*  w 


Jt=  1  n=l 
K  N 


=  *  -  \  E  E  li  IIJ,., 


(43) 


Jb=l  n=l 


where  6^  \tn)  is  the  conditional  mean  estimate  of  6{tn )  computed  at  9  =  9^l\  and  r^t;  •)  is 
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defined  by  (34)  -  (36). 


Thus,  to  compute  Q(0,0^)  (E -step)  we  only  need  to  compute  6^(t„)  n  =  1,2, . . .,  and 
the  maximization  of  Q(0, 6^)  (M-Step)  is  equivalent  to  the  minimization  of  each  of  the  terms 
in  the  k  sum  separately.  Therefore,  the  EM  algorithm  assumes  the  form: 

E-step:  Compute 

3(<,(tn)  =  Ed(I){^(tn)/z}  n—  1,2, ...  IV  (44) 

M-step:  For  k  =  1, 2, ...  K 


Min  £  ||  •(*)  -  lljw  — 

n=l  n 

Ok 


(45) 


The  algorithm  is  illustrated  in  Figure  5.  Once  again  we  note  that  the  minimization  in 
(45)  finds  the  ML  location  estimate  of  the  ktfi  signal  source  with  respect  to  the  estimated 
sensor  locations.  Thus,  the  algorithm  iterates  back  and  forth,  using  the  current  estimate  6^ 
of  the  various  source  location  parameters  to  improve  the  estimate  of  the  array  uncertainties 
and  thus  to  improve  the  next  source  location  estimates.  Since  the  algorithm  is  based  on  the 
EM  method,  it  is  ensured  to  converge  monotonically  to  the  joint  ML  estimate  of  all  source 
location  parameters,  or  at  least  to  a  stationary  point  of  the  joint  log-likelihood  function. 

All  the  effect  of  array  calibration  enters  the  E-step,  where  we  use  the  data  available  from 
all  signal  sources  to  form  the  best  possible  estimate  of  the  array  uncertainties.  We  note  that 
if  some  of  the  signal  sources  have  known  locations,  they  are  still  used  in  the  E-step  (with 
their  exact  locations).  In  the  M-step,  of  course,  we  only  re-estimate  the  unknown  location 
parameters. 

Perhaps  the  most  striking  feature  of  the  algorithm  is  that  it  decomposes  the  multi- 
dimensioned  optimization  associated  with  the  direct  ML  approach  into  optimizations  with 
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respect  to  each  source  location  parameters  separately  leading  to  a  considerable  simplification 
in  estimator  structure  and  computations.  Because  of  its  parallel  structure,  the  algorithm 


incorporates  any  pre-specified  number  of  signal  sources  with  only  a  linear  increase  in  the 
computational  complexity. 

To  compute  6^l\tn),  we  need  to  perform  the  conditional  expectation  indicated  in  (44). 
For  that  purpose,  we  concatenate  the  various  measurement  equations  in  (40): 


=  Hn  r(tn;0)  +  vn 

=*  Hnr{tn]0)--  HnDT{tn-  0)6(tn)  + 
c 

where  we  define: 

Zn  = 


r(f;0i) 

r(i;02) 

r{t]9K) 

r{t\9x) 

2) 

r{fJK) 
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(54) 


zn- Hnr(tn]e{l)) 

where  zn,  r  (t;  0),  D(t-,$),  Hn,  and  R„  are  given  by  (47),  (49),  (50),  (51),  and  (53),  respec¬ 
tively.  Since  6^\tn)  depends  only  on  the  measurements  zn  at  t  =  tn,  they  can  be  computed 
recursively  in  n.  There  is  no  need  to  store  past  measurements  for  the  purpose  of  estimating 
the  array  uncertainties. 
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V.  Gradient-Based  Algorithms  and  Performance  Evaluation 

As  indicated  in  (9j  (llj,  the  rate  of  convergence  of  the  EM  algorithm  is  exponential  (linear). 
Therefore,  near  the  point  of  convergence,  we  may  want  to  use  faster  converging  gradient-based 
methods  such  as  the  Newton-Raphson  or  the  Scoring  algorithm.  These  methods  require  the 
computation  of  the  log-likelihood  gradient,  and  the  computation  of  the  log-likelihood  Hessian 
or  the  Fisher’s  information  matrix  (FIM),  or  some  approximation  of  which.  The  FIM  can 
further  be  used  to  assess  the  m.s.e.  of  the  resulting  ML  parameter  estimates. 

The  computation  of  the  log-likelihood  gradient  (score)  dL%(0)/d0,  obtained  by  differ¬ 
entiating  (17),  is  extremely  complicated  since  it  involves  the  differentiation  of  the  Kalman 
filtering  equations.  An  alternative  approach  that  yields  the  same  result  with  much  simpler 
computations  is  based  on  Fisher’s  identity  [12]: 

Ail(9)  =  i;<(^iy(9)/z]  (55) 

Fisher’s  identity  asserts  that  the  observed  data  score  jgL%(9)  can  be  computed  by  taking  the 
conditional  expectation  of  the  complete  data  score  jqL y(0)  .  A  sketch  proof  of  the  identity 
is  presented  in  Appendix  B.  A  formal  proof  with  the  associated  regularity  conditions  can  be 
found  in  [13], 

Consider  first  the  single  source  case.  Differentiating  (28) 

=  f„r(e„;9))  (56) 

n=  1 

where  the  columns  dr£(t;  0)/50  of  drT{t\9)/d9  are  closely  approximated  using  Taylor 
series  expansion  about  the  nominal  sensor  locations: 


dr£(t;0) 

d9 

Differentiating  (2): 


_  dr£(t;g) 


dTm{t\0)  _  1  dpT(t;9) 
99  |q„(4)  c  d9 


<M«;*) 


(57) 

(58) 
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dpTd{oe)  *).<£(*;*)]  (59) 


Differentiating  (4) 

d2T*(t-,0)  _  1 _ 1 

c||P(t;0)-  hm  (t)  ||  ae 

where  d m(t;0)  is  the  unit  vector  defined  in  (9).  Substituting  (58)  and  (59)  into  (57)  and 

O 

ignoring  terms  of  the  order  of  ||  6m(t)  ||  /  ||  p(t;0)-  qm  (t)  ||  (that  is,  the  ratio  of  array 
uncertainties  to  the  source  range), 


d»m(M)  ldprM),  ,  -v  dTm{t-0) 

—do~  =  c~a s) - To — 


(60) 


Eq.  (60)  asserts  that  to  a  first  order  approximation,  the  sensor  displacements  have  negligible 
effect  on  the  partial  derivatives  drm(t;  0)/d0.  Therefore,  using  (60)  in  (56) 

yol^0)  -  £  tz*  - 


n=l 


n=l 


30 


zn  -  Hn  t  (t„;0)  +  -ffnZ?T(tn;0)*(tn) 


(61) 


Substituting  (61)  into  (55),  the  observed  data  score  is  given  by: 


life)  =  f;  ■?<>=*>  gjg-1 


n=l 


30 


zn-HnT  (tn;0)  +  -tf„DT(tn;0)<5(in) 
c 


(62) 


where  £(tn)  is  the  conditional  mean  estimate  of  0(fn): 


6{tn)  =  Ee{6{tn)/z)  (63) 

The  computation  of  tf(tn)  can  be  carried  out  efficiently  using  the  smoothing  equations 
presented  in  Appendix  A.  We  have  therefore  substituted  the  direct  differentiation  of  (17)  that 
involves  the  differentiation  of  the  Kalman  filtering  equations  with  a  much  simpler  computation 
that  only  involves  the  Kalman  smoothing  equations. 
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Taking  the  second  moment  (con variance)  of  the  score,  the  FIM  is  given  by: 


m 


Eg 


dLl{6)dL,{0)\ 
dO  dO  j 


N  M 

EE 

n=l  m=l 


of  ( tni0 ) 

dO 


HlR-'Eg  {v„v£}  R^Hr 


dr(tm\6) 
1  dO 


(64) 


where: 

vn  =  z n-Hnr  (t„;  0 )  +  - HnDT(tn ;  <?)$(£„)  (65) 

c 

Recursive  formulae  for  calculation  Eg  {vr»Vjt  j  can  found  in  [14], 

To  compute  the  log-likelihood  Hessian  we  must  differentiate  the  expression  for  the  score. 
In  practice,  the  Hessian  may  be  approximated  by  replacing  derivatives  with  finite  differences. 
To  approximate  the  partial  derivatives  of  the  score,  we  perturb  one  coordinate  of  0  at  a 
time  and  compute  the  resulting  score  according  to  (62).  This  approximation  requires  the 
computation  of  the  Kalman  smoothing  equations  at  closely  spaced  values  of  0,  a  task  that 
cam  be  simplified  by  pre-computation  of  the  smoothed  error  covariance  matrices.  Another 
approximation  of  the  Hessian  and  the  FIM,  based  on  score  computation,  is  presented  in  [14]. 

In  the  multiple  source  case,  the  score  and  the  FIM  are  still  given  by  (62)  and  (64),  where 
z„,  t  ( t;0 ),  D(t\0),  Hn,  and  Rn  are  given  by  (47),  (49),  (50),  (51),  and  (53),  respectively 
and  d  r  (t\0)/d0  is  the  corresponding  block  diagonal  matrix. 

The  formule  for  the  score  and  the  FIM  are  simplified  if  the  uncertainty  process  6{tn)  n  = 
1,2, .. .  is  uncorrelated  in  time.  Recall  Eq.  (39),  we  can  write  in  this  case 


^(£n) 


Eg{6(tn)/zn} 

-- QnD(tn-,0)H l  \±HnDT{tn]0)QnD{tn-0)Hl  +  R„ 
c  Lc* 

\zn  —  Hn  T  (<„;*)] 


(66) 


Substituting  (66)  into  (62)  and  following  straightforward  matrix  manipulations,  we  obtain 


the  following  explicit  formula  for  the  score: 


N  o*  .  _  ^ 

=  £  d--d-~~-Hn  ±HnDT(tn-,0)QnD{tn]0)HZ  +  Rn  [zn  -  Hn  r  (t„;*)] 

v  n=  1  L 


By  the  virtue  of  (15), 


zn  -  Hn  t  ( tn\0 )  ~  --HnDT{tn-0)6(tn)  +vn 


Substituting  (55)  into  (54)  and  taking  the  second  moment,  we  obtain  the  following  explicit 

formula  for  the  FIM: 

T 

AW  =  E  dT-^]  Hi  \±HnDT(tn]O)QnD(tn-,0)HZ  +  *«]  1  Hn— (69) 

Once  again  we  note  that  the  expression  in  (67)  and  (69)  are  valid  also  for  the  multiple 
source  case.  The  result  in  (69)  is  a  generalization  of  the  results  developed  in  [1],  [2],  [4], 
[5].  This  expression  can  be  used  to  analyze  the  attainable  m.s.e.  for  any  pre-specified  array- 
source  geometry  and  any  given  delay  measurement  scheme,  and  the  realizable  improvement 
in  localization  accuracy  obtainable  due  to  the  array  calibration  effect. 
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VI.  Conclusion 

We  have  developed  a  computationally  efficient  iterative  algorithm  for  source  localization 
and  tracking  using  active/passive  arrays  with  uncertainties  in  sensor  locations.  We  suppose 
that  the  available  data  consist  of  time  delay,  or  differential  time  delay,  measurements  of 
the  signal  wavefront  across  the  array.  We  consider  a  general  senario  in  which  the  array 
uncertainties  may  be  correlated  in  time  and  in  space.  The  proposed  algorithm  is  optimal 
in  the  sense  that  it  converges  montonically  to  the  Maximum  Likelihood  (ML)  estimate  of 
the  source  trajectory  parameters.  In  the  case  of  multiple  sources,  the  algorithm  makes  an 
essential  use  of  the  information  available  from  all  sources  to  reduce  the  array  uncertainties 
(tne  so-called  array  calibration)  and  thus  to  improve  the  localization  accuracy  of  each  signal 
source.  We  also  derived  new  expressions  for  the  log-likelihood  gradient,  the  Hessian,  and  the 
Fisher’s  information  matrix,  that  may  be  used  for  efficient  implementation  of  gradient  based 
algorithms,  and  for  assessing  the  mean  square  error  of  the  resulting  ML  parameter  estimates. 
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Figure  2:  Effect  of  Sensor  Uncertainties  on  Array-Source  Geometry. 
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Figure 


Figure  4 


Figure  4:  Array-Source  Geometry  for  Multiple  Sources. 
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Figure  5:  The  Proposed  Algorithm  for  Multiple  Sources. 
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Figure 


Appendix  A:  Kalman  Smoothing  Equations  for  Array 
Uncertainty  Estimation 


The  array  uncertainties  £(£„)  n  =  1,2,...  and  the  observed  data  z„  n  =  1,2, 
following  stochastic  dynamic  linear  state  equations: 

6(tn)  =  $n<5(*n-l)  +  U„ 

z„  =  Hn6(tn)+\n  n  =  1,2, ...  iV 

Define: 


#V*  =  Ee{6(tn)/z1z2---Zk > 

pn/k  =  Eg  {[«(*„)  [«(«n)-A*„/t]r/Zl,Z2-"ZJt| 

Propagation  Equations: 

Forn  =  1,2, ...  iV 

Pn/n-l  ^nPn—l/n—l 

Pn/n—  1  =  r»—  1/n-l^n  b  Qn 

Up-Dating  Equations: 

Forn  =  1,2,...  iV 

Mr»/n  =  Pn/n-l  b  Kn  —  HnHn/n_  jJ 
Pn/n  =  [I  ~  EnHn\  Pn/n-i 


I*  0/0  —  0 
Po/o  =  Po 


. . .  satisfy  the 
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where  Kn  is  the  Kalman  Gain: 


Kn  —  Pn/n-l^n  [-^n Pn/n-1 +  Pn\ 

Smoothing  Equations: 

Forn  =  N,N  - 

A*n— l/f»  /*n-l/n-l  "l"  &n—l  [l*n/N  ^nPn-l/n-l] 

Pn-l/N  ~  Pn-l/n-1  +  •S'n-l  \Pn/N  ~  ^n/n-l]  ^n-l 

where 

$n- 1  =  Prx-l/n-l^n^n/n-l 

The  outcome  of  the  smoothing  equations  are  the  uncertainty  estimates  and  the  associated 
error  covariances: 


*(*«)  =  Pn/N 

P{pn)  ~  P ifN 
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Appendix  B:  Fisher’s  Identity 


Let  i  and  y  =  ( zT:6T)T  denote  the  observed  (incomplete)  and  the  complete  data  vectors, 
respectively.  Using  Baye’s  rule  for  probability  densities, 

/y(y;0)  =  /z(z;0)  •  /a/z(^/z;0) 

Taking  the  logarithm  on  both  sides  of  the  equation  and  differentiating  with  respect  to  9 , 

^i°g/z(z;*)  =  ^i°g/Y(y;*)  -  -^i°g/A/z(Vz;*) 

Taking  the  conditional  expectation  given  z, 

j^l°g/z(z;0)  =  Ee  |^log/Y(y;tf)/zJ  -  E0  j  j^log/A/z(Vz;0)/z} 

Now, 

^{^log/A/z(Vz;«)/z}  = 

=  J^ogfA/z(6/z;8)-fA,z(6/z;0)dS  =  J  ^fA/z(S/z;0)dS 

=  Tel  f A, z(.6W)d*  =  §-eW  =  o 

Therefore, 

^i°g/z(z;fl)  =  j^iog/Y(y;0)/z} 

or, 

which  is  Fisher’s  identity. 


39 


DOCUMENT  LIBRARY 

July  5.  1989 

Distribution  List  for  Technical  Report  Exchange 


Attn:  Stella  Sanchez-Wade 
Documents  Section 
Scripps  Institution  of  Oceanography 
Library,  Mail  Code  C-075C 
La  Jolla.  CA  92093 

Hancock  Library  of  Biology  & 
Oceanography 
Alan  Hancock  Laboratory 
University  of  Southern  California 
University  Park 
Los  Angeles,  CA  90089-0371 

Gifts  &  Exchanges 
Library 

Bedford  Institute  of  Oceanography 
P.O.  Box  1006 

Dartmouth.  NS.  B2Y  4A2.  CANADA 

Office  of  the  International 
Ice  Patrol 

c/o  Coast  Guard  R  &  D  Center 
Avery  Point 
Groton,  CT  06340 

Library 

Physical  Oceanographic  Laboratory 
Nova  University 
8000  N.  Ocean  Drive 
Dania.  FL  33304 

NOAA/NESDIS  Miami  Library  Center 
4301  Rickenbacker  Causeway 
Miami.  FL  33149 

Library 

Skidaway  Institute  of  Oceanography 
P.O.  Box  13687 
Savannah,  GA  31416 

Institute  of  Geophysics 
University  of  Hawaii 
Library  Room  252 
2525  Correa  Road 
Honolulu.  HI  96822 

Library 

Chesapeake  Bay  Institute 
4800  Atwell  Road 
Shady  Side.  MD  20876 

MIT  Libraries 

Serial  Journal  Room  14E-210 
Cambridge.  MA  02139 


Director,  Ralph  M.  Parsons  Laboratory 

Room  48-311 

MIT 

Cambridge,  MA  02139 

Marine  Resources  Information  Center 

Building  E38-320 

MIT 

Cambridge,  MA  02139 
Library 

Lamont- Doherty  Geological 
Observatory 
Colombia  University 
Palisades,  NY  10964 

Library 

Serials  Department 
Oregon  State  University 
Corvallis,  OR  97331 

Pell  Marine  Science  Library 
University  of  Rhode  Island 
Narragansett  Bay  Campus 
Narragansett,  RI  02882 

Working  Collection 
Texas  A&M  University 
Dept,  of  Oceanography 
College  Station,  TX  77843 

Library 

Virginia  Institute  of  Marine  Science 
Gloucester  Point,  VA  23062 

Fisheries-Oceanography  Library 
151  Oceanography  Teaching  Bldg. 
University  of  Washington 
Seattle,  WA  98195 

Library 

R.S.M.A.S. 

University  of  Miami 

4600  Rickenbacker  Causeway 

Miami,  FL  33149 

Maury  Oceanographic  Library 
Naval  Oceanographic  Office 
Stennis  Space  Center 
NSTL,  MS  39522-5001 

Marine  Sciences  Collection 
Mayaguez  Campus  Library 
University  of  Puerto  Rico 
Mayagues,  Puerto  Rico  00708 


Mac  88- 11 6  (*1 6B) 


50272-101 


REPORT  DOCUMENTATION  1  REPORT  NO.  2. 

PAGE  WHOI-89-32 

3.  Recipient's  Accession  No. 

4.  THIe  and  Subtitle 

Optimal  Source  Localization  and  Tracking  Using  Arrays  with  Uncertainties  in  Sensor 
Locations 

1  Report  Deis 

August  1989 

6. 

7.  Authors) 

Mordechai  Segal  <nid  Chad  Weinstein 

8.  Performing  Organization  Rept.  No. 

WHOI-89-32 

9.  Performing  Organization  Name  and  Addraes 

The  Woods  Hole  Oceanographic  Institution 

Woods  Hole,  Massachusetts  02543 

10.  Project/Taak/Worii  Unit  No. 

11.  Contract(C)  or  Qrant(G)  No. 

(C)  N00014-85-K-0272 
«J) 

12.  Sponaoring  Organization  Name  and  Addraaa 

The  Office  of  Naval  Research 

13.  Type  of  Report  A  Period  Covered 

Technical  Report 

14. 

15.  Supplementary  Notes 

This  report  should  be  cited  as:  Woods  Hole  Oceanog.  Inst.  Tech.  Rept.,  WHOI-89-32. 

16.  Abstract  (Limit:  200  words) 

We  develop  a  computationally  efficient  iterative  algorithm  for  source  localization  and  tracking  using  active/passive  arrays  with 
uncertainties  in  sensor  locations.  We  suppose  that  the  available  data  consist  of  time  delay,  or  differential  time  delay,  measurements  of 
the  signal  wavefront  across  the  array.  We  consider  a  general  scenario  in  which  the  array  uncertainties  may  be  correlated  in  time  and  in 
space.  The  proposed  algorithm  is  optimal  in  the  sense  that  it  converges  montonically  to  the  Maximum  Likelihood  (ML)  estimate  of  the 
source  trajectory  parameters.  In  the  case  of  multiple  sources,  the  algorithm  makes  an  essential  use  of  the  information  available  from  all 
sources  to  reduce  the  array  uncertainties  (the  so-called  array  calibration)  and  thus  to  improve  the  localization  accuracy  of  each  signal 
source.  We  also  derive  new  expressions  for  the  log-likelihood  gradient,  the  Hessian,  and  the  Fisher's  information  matrix,  that  may  be 
used  for  efficient  implementation  of  gradient  based  algorithms,  and  for  assessing  the  mean  square  error  of  the  resulting  ML  parameter 
estimates. 

17.  Document  Analysis  a.  Descriptors 

1.  source  localization  4.  EM  Algorithm 

2.  sensor  uncertainties  5.  Array  calibration 

3.  Maximum  Likelihood 

b.  IdentlflerWOpen-Ended  Terms 


e.  COSAT1  ReWOroup 


18.  Availability  Statement 

19.  Security  Class  (This  Report) 

21.  No.  of  Pages 

Approved  for  publication:  distribution  unlimited. 

UNCLASSIFIED 

39 

20.  Security  Class  (This  Pegs) 

22.  Price 

(See  AN8I-Z39.18)  See  In tt met  loot  on  Rtvru  OPTIONAL  FORM  272  (4-77) 

(Formerly  NTIS-35) 
Department  of  Commerce 


